#include "rasterToField.hpp"
#include "fieldFunctions.hpp"

namespace hms::raster
{

scalar rasterValueAtCoords(
	const Raster& raster,
	const Ref<const Array2s>& coords
){
	/* cast to int instead of floor() is ok, because 
	 * coords -  raster.bottomLeft() is always positive when in bounds, 
	 * and out of bounds indices are caught afterwards */
	Array2i index {
		( ( coords - raster.bottomLeft() ) / raster.cellsize() ).cast<Index>()
	};

	/* make sure indices are in bounds */
	index = (index < 0).select(0, index);
	index.x() = std::min( index.x(), raster.nx()-1 );
	index.y() = std::min( index.y(), raster.ny()-1 );

	return raster.values()( index.x(), index.y() );
}

template<typename MeshType>
ArrayXXs rasterToFieldValues(
	const Raster& raster,
	const MeshBase<MeshType>& mesh
){
	ArrayXXs fieldVals;
	resizeFieldValues( fieldVals, 1, mesh, NonUniform );

	mesh.forEachCell( [&](const auto& cell){
		fieldVals( 0, cell.index() ) = rasterValueAtCoords(
			raster, cell.centroid()
		);
	});

	return fieldVals;
}

#define INSTANTIATE_FOR_ALL_MESHTYPES(MeshType) \
template ArrayXXs rasterToFieldValues( const Raster&, const MeshBase<MeshType>& );
#include "meshInstantiator.hpp"
	
} // namespace hms::raster
